# Load Data ---------------------------------------------------------------
# Set file path
setwd("")

load("drop_in_ocean_data.rdata")

# Load Required hetreg Function --------------------------------------------------
source("09_hetReg.R")


# Analyses ----------------------------------------------------------------
# * By Race ---------------------------------------------------------------
# estimate model among Black respodents
m_cjs_att_blk <- lm(cjs_att ~ direct_exp + rr_sc + 
                      peer.felony + age_sc + man + 
                      educ_sc + inc_sc + pid7_dem,
                    data = cjs.df, weights = wts_bw_stk, black == 1)
summary(m_cjs_att_blk)

# store data frame from models
dat_b <- as.matrix(m_cjs_att_blk$model)

# estimate heteroskedastic regression
m_cjs_att_blk_var <- MLhetreg(dat_b[,c("cjs_att")], 
                              dat_b[,c(2:(ncol(dat_b)-1))], 
                              dat_b[,c("direct_exp","rr_sc", "peer.felony",
                                       "age_sc", "man", "educ_sc", "inc_sc", "pid7_dem")]) 
summary.MLhetreg(m_cjs_att_blk_var)


# estimate model among White respodents
m_cjs_att_wht <- lm(cjs_att ~ direct_exp + rr_sc + 
                      peer.felony + age_sc + man + 
                      educ_sc + inc_sc + pid7_dem,
                    data = cjs.df, weights = wts_bw_stk, black == 0)
summary(m_cjs_att_wht)

# store data frame from models
dat_w <- as.matrix(m_cjs_att_wht$model)

# estimate heteroskedastic regression
m_cjs_att_wht_var <- MLhetreg(dat_w[,c("cjs_att")], 
                              dat_w[,c(2:(ncol(dat_w)-1))], 
                              dat_w[,c("direct_exp","rr_sc", "peer.felony",
                                       "age_sc", "man", "educ_sc", "inc_sc", "pid7_dem")]) 
summary.MLhetreg(m_cjs_att_wht_var)


# * By Gender (Black) -----------------------------------------------------
# estimate model among women respodents
m_cjs_att_blk_wom <- lm(cjs_att ~ direct_exp + rr_sc + 
                          peer.felony + age_sc + educ_sc + inc_sc + pid7_dem,
                        data = cjs.df, subset = black == 1 & man == 0, weights = wts_bw_stk)
summary(m_cjs_att_blk_wom)

# store data frame from models
dat_bw <- as.matrix(m_cjs_att_blk_wom$model)

# estimate heteroskedastic regression
m_cjs_att_blk_wom_var <- MLhetreg(dat_bw[,c("cjs_att")], 
                                  dat_bw[,c(2:(ncol(dat_bw)-1))], 
                                  dat_bw[,c("direct_exp","rr_sc", "peer.felony",
                                            "age_sc",  "educ_sc", "inc_sc", "pid7_dem")]) 
summary.MLhetreg(m_cjs_att_blk_wom_var)


# estimate model among men respondents
m_cjs_att_blk_man <- lm(cjs_att ~ direct_exp + rr_sc + 
                          peer.felony + age_sc + educ_sc + inc_sc + pid7_dem,
                        data = cjs.df, subset = black == 1 & man == 1, weights = wts_bw_stk)
summary(m_cjs_att_blk_man)


# store data frame from models
dat_bm <- as.matrix(m_cjs_att_blk_man$model)

# estimate heteroskedastic regression
m_cjs_att_blk_man_var <- MLhetreg(dat_bm[,c("cjs_att")], 
                                  dat_bm[,c(2:(ncol(dat_bm)-1))], 
                                  dat_bm[,c("direct_exp","rr_sc", "peer.felony",
                                            "age_sc", "educ_sc", "inc_sc", "pid7_dem")]) 
summary.MLhetreg(m_cjs_att_blk_man_var)

# * By Gender (White) -----------------------------------------------------
# estimate model among women respodents
m_cjs_att_wht_wom <- lm(cjs_att ~ direct_exp + rr_sc + 
                          peer.felony + age_sc + educ_sc + inc_sc + pid7_dem,
                        data = cjs.df, subset = black == 0 & man == 0, weights = wts_bw_stk)
summary(m_cjs_att_wht_wom)


# store data frame from models
dat_ww <- as.matrix(m_cjs_att_wht_wom$model)
# estimate heteroskedastic regression
m_cjs_att_wht_wom_var <- MLhetreg(dat_ww[,c("cjs_att")], 
                                  dat_ww[,c(2:(ncol(dat_ww)-1))], 
                                  dat_ww[,c("direct_exp","rr_sc", "peer.felony",
                                            "age_sc",  "educ_sc", "inc_sc", "pid7_dem")]) 
summary.MLhetreg(m_cjs_att_wht_wom_var)


# estimate model among men respondents
m_cjs_att_wht_man <- lm(cjs_att ~ direct_exp + rr_sc + 
                          peer.felony + age_sc + educ_sc + inc_sc + pid7_dem,
                        data = cjs.df, subset = black == 0 & man == 1, weights = wts_bw_stk)
summary(m_cjs_att_wht_man)

# store data frame from models
dat_wm <- as.matrix(m_cjs_att_wht_man$model)
# estimate heteroskedastic regression
m_cjs_att_wht_man_var <- MLhetreg(dat_wm[,c("cjs_att")], 
                                  dat_wm[,c(2:(ncol(dat_wm)-1))], 
                                  dat_wm[,c("direct_exp","rr_sc", "peer.felony",
                                            "age_sc",  "educ_sc", "inc_sc", "pid7_dem")]) 
summary.MLhetreg(m_cjs_att_wht_man_var)


